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and Kate Tilling b 

Growth models are commonly used in life course epidemiology to describe growth trajectories and their deter- 
minants or to relate particular patterns of change to later health outcomes. However, methods to analyse 
relationships between two or more change processes occurring in parallel, in particular to assess evidence for 
causal influences of change in one variable on subsequent changes in another, are less developed. We discuss 
linear spline multilevel models with a multivariate response and show how these can be used to relate rates 
of change in a particular time period in one variable to later rates of change in another variable by using the 
variances and covariances of individual-level random effects for each of the splines. We describe how regression 
coefficients can be calculated for these associations and how these can be adjusted for other parameters such 
as random effect variables relating to baseline values or rates of change in earlier time periods, and compare 
different methods for calculating the standard errors of these regression coefficients. We also show that these 
models can equivalently be fitted in the structural equation modelling framework and apply each method to 
weight and mean arterial pressure changes during pregnancy, obtaining similar results for multilevel and struc- 
tural equation models. This method improves on the multivariate linear growth models, which have been used 
previously to model parallel processes because it enables nonlinear patterns of change to be modelled and the 
temporal sequence of multivariate changes to be determined, with adjustment for change in earlier time periods. 
Copyright © 2012 John Wiley & Sons, Ltd. 
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1. Introduction 

Growth models are often used in life course epidemiology to describe patterns of growth or change, 
to assess the influences of rates and timing of change on later health outcomes and to investigate the 
effects of exposures on patterns of change [1-6]. In some cases two or more change processes occur in 
parallel, such as gains in height and weight during childhood, and it may be of interest to assess rela- 
tionships between such processes, or to adjust for the effects of one change process when assessing the 
associations of another with a later health outcome. It may also be hypothesised that change in one pro- 
cess causes a change in another, or that an unmeasured or latent construct causes change in both of the 
processes simultaneously. Bivariate multilevel models (MLMs) for repeated measurements with linear 
relationships of both outcome variables with time have been used to demonstrate the correlation between 
values of the two variables at baseline and rates of change in the variables over time [7-9]. However, this 
does not demonstrate whether change in one variable temporally precedes change in another and is thus 
unable to provide evidence towards a causal effect. 
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Linear spline models have recently been proposed as a method of representing growth [5] or change 
[10], which reduces the dimensionality of repeated measurements. The shape of the trajectory of change 
is assumed to be piecewise linear, with knot points defining changes in the magnitude or direction of 
association of the response variable with time. The selection of the number and location of knot points 
may be determined by the data or by prior knowledge and therefore the linear splines represent the shape 
of the change trajectory in a meaningful way, and the linear periods of change provide coefficients that 
are more interpretable than for polynomial curves. By first fitting splines to each change process and 
then modelling these processes in parallel, we are able to model associations between observed rates 
of change in one variable in one period of time and observed rates of change in another variable in a 
subsequent period of time, with time periods defined by the splines, and thus we may identify whether 
changes in one variable precede changes in another. By examining the temporal sequence of events 
there is the possibility of improving causal inferences for pathways between variables. Furthermore, 
prior knowledge or hypotheses about which pathways exist may be tested by constraining covariances 
between individual-level random effects for different periods of change to be zero and comparing the 
model fit to models where these covariances are freely estimated. 

We will use the example of weight and mean arterial pressure (MAP) during pregnancy to demonstrate 
these models. In normal pregnancy there is a decrease in blood pressure in early pregnancy followed by 
a rise in late pregnancy [11]. Hypertensive disorders of pregnancy (HDP), defined by high blood pres- 
sure in late pregnancy (after 20 weeks' gestation), are associated with risk of adverse health outcomes 
for both the mother and offspring [12-14]. and gestational weight gain (GWG) has been found to be 
positively associated with the risk of developing an HDP [15-17]. However, it is not clear whether 
changes in weight during pregnancy precede changes in blood pressure, or whether increases in blood 
pressure precede GWG because of increases in oedema and plasma volume expansion [18]. MAP is 
the average pressure in an artery over a complete cycle of one heart beat. It is estimated by combining 
systolic and diastolic blood pressure and allowing for the lower pressure during the diastolic phase of 
the cardiac cycle (see below for description of how this is calculated). In this paper we have used MAP 
so that we have one blood pressure variable that we can relate to weight change in pregnancy and that 
takes account of blood pressure in both systole and diastole. We show how MLMs and structural equa- 
tion models (SEMs) with linear splines for gestational age can be applied to weight and MAP to learn 
whether an increase in weight precedes a rise in MAP, and in which periods of pregnancy associations 
are strongest. The paper is structured as follows: in Section 2 we describe the data. In Section 3 we 
define the multilevel multivariate response model and the derivation of regression coefficients from the 
variances and covariances of random effects. In Section 4 we show that the model can be equivalently 
defined in SEM form. In Section 5 we provide an application of the model to weight and MAP changes in 
pregnancy, and compare results from fitting the model in MLM software (MLwiN) (Rasbash J, Browne 
W, Healy M, Cameron B, Charlton C, Centre for Multilevel Modelling, University of Bristol, UK) and 
SEM software (Mplus) (Muthen and Muthen, Los Angeles, California). Section 6 is a discussion of the 
methods described. 



2. Data 



The data are from the Avon Longitudinal Study of Parents and Children (ALSPAC), which is described in 
full elsewhere [19] and on the website http://www.bris.ac.uk/alspac. A total of 14,541 pregnant women 
were recruited, who were living in a defined area of Avon including the city of Bristol during their preg- 
nancy and had an expected delivery date between 1 April 1991 and 31 December 1992. Ethical approval 
for the study was obtained from the ALSPAC Law and Ethics Committee and from the National Health 
Service local ethics committee. We restrict analysis to singleton live term births (^ 37 weeks' gesta- 
tion) with no evidence of pre-eclampsia (defined according to the International Society for the Study 
of Hypertension in Pregnancy criteria [20]) or a previous diagnosis of hypertension (N — 11, 650). All 
weight and blood pressure measurements, which were taken routinely as part of antenatal care by mid- 
wives or obstetricians, were abstracted from obstetric records by six trained research midwives. This 
resulted in a median of 14 (interquartile range 1 1 to 16) blood pressure measurements and median of 12 
(interquartile range 10 to 14) weight measurements per woman. MAP was calculated as 1/3 x systolic 
blood pressure + 2/3 x diastolic blood pressure. Confounding variables considered here are maternal 
height, age, parity, education, smoking and offspring sex. Maternal age and offspring sex were obtained 
from obstetric records and maternal height, parity, highest educational qualification and smoking sta- 
tus were obtained from questionnaires administered during pregnancy. Smoking status was classified as 
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Table I. Characteristics of the dataset (N = 11, 650) and the subset with complete data on all covariates 
(N = 9429). 


Maternal characteristic 


A' with 
data 


Mean (SD) or % 


Mean (SD) or % in 
complete case 


Height (cm) 


10,278 


164.0 (6.69) 


164.1 (6.65) 


Age (yrs) 


11,650 






<20 




4.84 


3.45 


20-24 




19.30 


16.66 


25-29 




38.66 


39.56 


30-34 




27.60 


29.95 


35+ 




9.60 


10.38 


Parity 


10,830 






Nulliparous 




44.56 


44.63 


Multiparous 




55.44 


55.37 


Smoking in pregnancy 


10,945 






Never 




66.20 


68.69 


Prepregnancy/first trimester 




13.85 


13.41 


Throughout 




19.95 


17.90 


Highest educational qualification 


10,492 






CSE/vocational 




29.89 


28.10 


0 level 




34.51 


34.99 


A level 




22.57 


23.35 


Degree 




13.03 


13.55 


Sex of offspring 


11,650 






Male 




51.12 


50.79 


Female 




48.88 


49.21 



'never', 'immediately prepregnancy or first trimester only', or 'throughout pregnancy' for women who 
continued to smoke after the first trimester. Table I shows the characteristics of the women included in 
the analysis. Linear spline random effects models have previously been fitted to weight measurements 
in pregnancy, to model changes in weight with gestational age [21]. A similar method was used to fit a 
linear spline random effects model to MAP (see Section 5). 

3. Multilevel spline model with a multivariate response 



3.1. Multivariate response linear spline random effects model 

Multilevel models for repeated measurements have been described by Goldstein [22] and others. We 
will first consider a univariate MLM for repeated measurements of a response variable y, measured on 

occasions, j — 1, . . . , r& for each individual k — 1, n, with two levels: measurement occasion j 

within individuals . The response variable is modelled as a linear function of time and the intercept and 
linear slopes are allowed to vary at the individual level by including individual-level random effects for 
these parameters as follows: 

yjk = Pok + (Pik + Uik)t]k + u ok + £ojk (1) 

where y ^ is the value of the response variable, y, at the j th measurement for the k th individual, tjk is 
the time of the j* measurement for the k lh individual, « 0/ t and are the individual level (level-2) 
residuals for the intercept and slope respectively and the £ojk s are the measurement occasion level 
(level- 1) residuals. These parameters have the following assumptions: 



Uok 
Ulk 



N(0, £2, 



2 

eO 



'ul 



(2) 



A bivariate response model to investigate linear change in two variables for i — 1,2 over time 
has three levels: responses^ 1 ' within measurement occasion j within individual^ and thus a model with 
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a random intercept and random slope at the individual level is denned by the following formula and 
assumptions [22]: 
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(3) 



To extend this to a model where time is represented by a set of linear splines for response (i) 
where, for splines there are + 1 time points, t$ < fj < . . . < t m then 



hjk - 
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Hence, a linear spline model with a bivariate response has the form 

„,«■> 

3$ = A?' + E (#° + "« ) '«* + <' + for / = 1,2 



,(i) 
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r (D(2) 



r 2(l) 



r 2(2) 



(m< | )+m< 2 >+2)x(m<l>+m< 2 >+2) 



(5) 

The individual-level random effects give each individual's deviance from the average intercept (iiok) 
and from the average slope of each of the splines (wit, . . ., w m (Ofc) and the individual-level covariance 
matrix contains the pair wise covariances between each set of these random effects for the intercept and 
slopes within each of the response variables and between the response variables. For simplicity we have 
defined a bivariate response model, but this can be extended to more than two response variables. In this 
case the covariance matrix at the individual level would include all variances and covariances between 
individual-level residuals and the covariance matrix at the occasion level would include all variances and 
covariances between occasion-level residuals. 

Of note, the model estimated variances and covariances of individual-level random effects, a 2 ,^, 

6j^/' 2 \ are n °t the same as the variances and pairwise covariances of the estimated individual-level 
residuals produced by MLMs, uf^, W/^> because the residuals are shrunken towards the population 
mean by varying degrees, reflecting the level of uncertainty with which each individual's residual is 
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estimated [23]. Because the individual-level variance-covariance matrix estimated by the model pro- 
vides estimates of the population variances and covariances of the random effects this should be used to 
assess associations between the intercepts and slopes in the population rather than the shrunken residuals. 

3.2. Testing hypotheses about temporal relationships between response variables 

The model defined in Equation (5) allows the individual-level random effects relating to the intercepts 
and each of the slopes of the response variables in different time periods to be correlated with each other. 
Specific hypotheses about time periods in which changes in each of the variables are associated may 
be applied by setting the covariances between the corresponding pairs of the random effects to be zero, 
thus simplifying the model. The hypotheses may be tested by comparing the simplified models with the 
full unrestricted model. The fit of the models may be compared by using a likelihood ratio test, with the 
difference in degrees of freedom between the models being the number of constrained covariances. 

For simplicity of explanation, in the following example hypotheses it will be assumed that the same 
splines are used for both response variables; however, these definitions can easily be adapted to take into 
account time periods with different endpoints. 

Hypothesis 1: Changes in each of the variables are only associated in the same or adjacent 
time periods 

This hypothesis can be tested by fitting a model with covariances between individual-level random 
effects relating to change in nonadjacent time periods to be zero: 

ce !hh > = 0 for all |/i -/ 2 | > 1; /i,/ 2 ^0, 

and comparing the fit of this model with the full model with no constraints. The hypothesis can be used 
to test whether change in one variable in one time period has long-term associations with change in the 
other variable, or whether it is only associated with change in the immediately subsequent time period. 

Hypothesis 2: Changes in variable 1 are associated with subsequent changes in variable 2, but 
changes in variable 2 are not associated with subsequent changes in variable 1 

For this hypothesis, in the restricted model, all covariances between random effects relating to change in 
variable 2 and random effects relating to change in variable 1 in a later time period are set to zero: 

°uhh ) = 0 for aU /^/a; h,l 2 ^0 

This can be used to test for evidence that changes in variable 1 precede changes in variable 2, rather 
than changes in variable 2 preceding change in variable 1. 

Hypothesis 3: There is a lag time between change in one variable and change in another 

To test this hypothesis, in the restricted model the covariance between random effects relating to change 
in one variable in a particular time period and random effects relating to change in the other variable in 
subsequent time periods is set to zero for any time period where it is hypothesised that there may be a 
lag, but covariances with random effects relating to change in time periods after the hypothesised lag 
time remain unconstrained 

a thh ) = Q fork-/i<S j; h,l 2 ^0, 

where S s is the number of splines/time periods that make up the hypothesised lag time. 

The choice of hypotheses to be tested should be guided by the research question, which the analysis is 
designed to answer, and prior knowledge about the plausible relationships between the variables studied. 

3.3. Deriving regression coefficients from the random effects variance-covariance matrix 

From the covariance matrix of the individual-level random effects relating to baseline values of the 
response variable or rates of change in these variables in each time period described by the splines it is 
possible to derive regression coefficients to describe the associations between baseline values and rates 
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:t M;' 1 ^ associated with a unit change in random effect uf^ is given by [24] 



of change in each of the variables in different time periods. The regression coefficient for the change in 

(i 
J, 

h,i 2 =l,2; Zi = 0,l,...,m (il) ; l 2 = 0, 1, . . .,m (h) (6) 
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This formula can be extended to give regression coefficients for the relationship between two variables 
when adjusting for others. Fisher described how to adjust regression coefficients for additional variables 
[25] and expressions for regression coefficients when adjusting for up to four covariates are given in the 
online supplemental material 1 . 

In the more general case, to regress random effect wj' 0 ^ on p independent random effect 



variables: u^uf^, . . ., uf p \ 
cients relating to uf^ u^ 2 \ . . 
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let /? = (fii, fo, ■ ■ ■, Pp)' be the p x 1 vector of regression coeffi- 

. . ., uf^\ respectively, adjusting for each of the other variables. Also, let 

(foWp) \ ' b e the n x 1 vector of covariances between and each of the 
ulo'p ) 'o 
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a. 
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a. 
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be the p x p matrix of variances and pairwise covariances of the independent variables. Then E/? = oq, 
meaning that the mutually adjusted regressions can be found by solving the equation: 



0 = E-Vo 

3.4. Standard errors of regression coefficients derived from the random effects 
variance-covariance matrix 



(7) 



We consider a number of methods for calculating the standard errors of the regression coefficients 
derived from the random effects variance-covariance matrix and compare the performance of these 
methods on the example dataset in Section 5. The first method estimates the standard errors as if we had 
sampled n observations of the individual-level random effects from an underlying population (where 
n is the number of individuals in our MLM), whereas Methods 2 and 3 estimate the standard errors 
by treating the regression coefficients as nonlinear combinations of the estimates of the variances and 
covariances of the random effects, but do this in different ways. 

Method 1: If we treat our estimates of the variances and covariances of the random effects as if they 
were obtained from a sample of the individual-level random effects, then the standard errors of the 
regression coefficients can also be derived using formulae described by Fisher [25]. For the general case 
with p independent variables, where /?, E and cto are as in (7) above, let d be the p x 1 vector containing 
the diagonal elements of the inverse of E, d — (E _1 [l, 1], E _1 [2, 2], . . ., E _1 [/?, p]) '. Then the p x 1 
vector of standard errors, S E p — (SE(/3i), SE(fj 2 ), . . ., SE{f5 p ))' is given by 



SE* = 



1 



(8) 



Method 2: The second approach is to use the delta method to approximate the standard errors [26, 27]. 
This can be implemented in Stata (StataCorp LP, College Station, Texas) using the command nlcom 
and by explicitly defining the formulae for the regression coefficients from the estimated variances and 
covariances of random effects (see supporting information for formulae)*. We note that some programs 



^Supporting information may be found in the online version of this article. 
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estimate variances and covariances on transformed scales [28], and implementations of the delta method 
are also available in R (The R Foundation for Statistical Computing, Vienna, Austria) [29]. This method 
is practical when adjusting for up to three or even four variables; however, because there is no convenient 
expression for the inverse of a square matrix of order greater than four, a method that solves the inverse 
of E numerically (see Method 3) is more appropriate when adjusting for more variables. 



Method 3: The third approach is to generate a large number, q,of realisations of a (m^ + + 
2) x 1 vector, which we shall call a, containing all of the variances and covariances of the ran- 
dom effects from the multivariate multilevel spline model (as described in Equation (5)), &h — 

. Ou0i h - °lih >■■■> ^ol|) o^co V for A = l,2,...,q, using a multivariate normal dis- 
tribution, a ~ N(y,F). The mean vector, v, of the multivariate normal distribution is the (m^ + 
m' 2 ' + 2) x 1 vector of estimates of the variances and covariances of the random effects from the 

MLM, v = (o%\ frffi,^, . . . *™) the («C0 + m« + 2) x (™<« + B « + 2) 

variance-covariance matrix, T, is the matrix of the variances and covariances of these estimates, 

var(ff£o 1} ) 



/»2(1) 

cov^q ,a u01 J var(a„ 01 J 



/~2(1) -2(2) \ /«2(2) \ 

After producing # realisations of a , the coefficients relating to the regression of one of the random 
effects on a subset p of the random effects from the MLM are derived for each of the realisations of a 
separately. The h th realisations of E, E/,, and <x 0 , croft, defined in (7) are formed from the estimates in a/, 
and the equation fi h — ET^tr q/, is solved for each of the q realisations. The overall estimate of the vector 

of regression coefficients /? = (pi,Pz> ■ ■ •> Pp^J is produced by taking the mean of the q realisations, 
fih = (f)i h , ^2 h , • ■ ■, Pp h j of each of these regression coefficients, ft — Ylt=i (fih/q^J- 95% confi- 
dence intervals for each of the elements of fi : /J 2 , ■ ■ ., fi p are then produced by using the 2.5 th and 
97.5 th percentiles of their distributions over the q realisations as the lower and upper limits, respectively. 

We have uploaded Stata commands that implement Methods 2 and 3, named ref f adj ust4nlcom 
and ref f adjustsim, respectively, to the Statistical Software Components archive [30]. 



4. Comparison with a structural equation model 

Structural equation models contain a measurement equation and a structural equation [31, 32]. In 
the measurement equation, each of the measurements yj loads on the latent variables fj where 
/ = 0, 1, . . .,m^\ The intercept term for variable i is and the loadings for each of the mea- 
surements yj' on this latent variable are 1. The other factors represent the slopes for each of the 

periods defined by the splines. Using the same + 1 time points, < tf^ < . . . <tm as above, the 
measurement equation is 



y? = ff> + E E -'&)//» 

htf^tj l:tfl 1 <tj<tf ) 
+ E °fi {i) + E f for; = 1,2. (9) 

Any rescaling of this equation proportional to these loadings can also be used. To see that this is 
equivalent to Equation (5), an extra subscript representing individuals, k, is added and \ &f^ + M /jt ) * s 
substituted for each of the ffP. This gives 
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yg = (/>S°+"£)+ E (f-fi)W 0+ ^) + E 
+ E o(^ +M {?) +£ « 

= + 4 } ) + E $ (#° + 4 } ) + £ & fo - = i. 2 - 
1=1 

(10) 

where sj^ is defined as in Equation (4). 

The structural equation of the SEM defines the mutual relationships between the factors. It can be 
shown that the regression coefficient for the univariate relationship between two factors, /, regressed 

' 1 

on is equivalent to that between two MLM random effects u^ 1 ^ and uf^ (given in (6)) by 

substituting (fif* + wf ) f° r each of the as follows: 



cov / f f) cov((^) + ^)),(^ + ^)) 

v*(f£°) var^f +«%>) 



Adjusted regression coefficients for the relationships between factors and their standard errors can be 
shown to be equivalent to those between random effects in a similar way. 

The MLM described in Equation (5) is equivalent to the SEM shown in the path diagram in Figure 1. 
There have been a number of papers demonstrating the equivalence of linear MLMs and SEMs when 
each individual has the same number of measurements occurring at the same time points meaning that 
the data are balanced [31-33]. It has also been shown that SEMs can be estimated equivalently to MLMs 
when data are unbalanced, by treating the design as balanced with missing data and using full informa- 
tion maximum likelihood estimation [33]. However, nonlinear MLMs and generalized linear MLMs do 
not always have equivalent parametrisations within the SEM framework [7,33]. For our weight and MAP 
example, where individuals have different numbers of measurements and these were taken at different 
times for each individual, this can be modelled in an SEM by including all times when measurements 
were taken on any individual (rounded to whole week periods) and treating a particular individual's 
measurements as missing for any time points on which they were not measured. 



5. Example: bivariate model of weight and MAP change in pregnancy 

Where there was more than one blood pressure or weight measurement in a 2-week period, one was 
selected at random for the analysis to prevent any individual from having too high an influence on the 
models. This left a median of 10 and range of 1 to 18 measurements per woman for analysis. Linear 
spline models were fitted to weight and MAP separately to describe the shapes of the average trajec- 
tories of weight and MAP with gestational age [11,21]. Full details of this process are provided in the 
online supporting information. The model that minimised the differences between predicted and actual 
values over pregnancy for weight had two knot points, at 18 and 28 weeks' gestation, and for MAP had 
three knot points, at 18, 30 and 36 weeks' gestation. For simplicity, because the first two knot point 
locations were similar for weight and MAP and altering the position of the second knot point by 1 week 
had little effect on the fit of these models, we have used knot points of 18 and 29 weeks for weight and 
18, 29 and 36 weeks for MAP. Baseline was set at 8 weeks' gestation for both weight and MAP because 
there were few measurements prior to this and model predictions earlier in pregnancy than this may 
be unreliable. Gestational age was rounded to the nearest whole week to make the MLM comparable 
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y[ 1] 




(2) 

y 2 







Figure 1. Path diagram for bivariate response spline model in structural equation format. Path diagram assumes 
that measurements y^P occur at equally spaced time points. 



with the SEM. The formulae for the separate multilevel spline models for weight and MAP, without any 
confounders included in models were 

Weight* = 62.943 + (o.314 + sfi* + (o.529 + + (o.461 + u*»*) 

weight weight 
+ U 0k "+" S Ojk 

MAP Jk = 81.572 + (-0.187 + < fc AP ) + (0.115 + t%£*) + (0.314 + uf k AP ) sf AP 



\jk t "Ok T" "Ojk ' 



(i.i82 +uT , )^+uir+e 0 r 



where s™ eight is the first gestational age spline for weight, up to 18 weeks' gestation, s^ 6 ' 811 ' is the second 
spline, from 18 — 29 weeks, and s™ eight is the third spline, from 29 weeks onwards. For MAP s^ AP is the 
first gestational age spline, up to 18 weeks' gestation, ^^ AP is the second spline, from 18 — 29 weeks, 
5^ AP is the third spline, from 29 — 36 weeks and s^ AP is the fourth spline, from 36 weeks onwards. 
The shape of the average weight and MAP trajectories across gestation predicted from these models are 
shown in Figure 2 , and Figure SI shows example predicted trajectories of weight and MAP for five 
individual women. 

We combined both of these separate spline models into one bivariate spline model, with two responses: 
weight and MAP. The distributions of weight and MAP at 8 weeks and rates of change in weight and 
MAP in each period of pregnancy predicted by this unadjusted bivariate model are shown in Table SI. 
We then included the potential confounders: maternal height (continuous variable); age (< 20, 20 — 24, 
25 — 29, 30 — 34, 35 + years); parity (nulliparous, multiparous); smoking (never, prepregnancy/first 
trimester, throughout pregnancy); education (CSE/vocational, O level, A level, degree) and offspring sex 
(male, female). Maternal height was centred about the mean (164 cm) and the reference categories for 
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Figure 2. Average trajectories of (a) weight and (b) mean arterial pressure across pregnancy predicted by 
univariate multilevel linear spline models (N = 11, 650). 



the other variables were: age 25 — 29 years, nulliparous, never smoked, O level qualification and male 
offspring. We completed analyses using Stata version 11.2 and MLwiN version 2.24 with runmlwin 
[34] for the MLMs and Mplus version 5 for the SEMs. 

The level 2 (individual-level) covariance matrix and correlations between individual-level random 
effects for this model, with all covariances between random effects freely estimated, are shown in 
Table II. We will call this MLM Model 1. We then considered several restricted models and tested 
different hypotheses about the relationships between changes in weight and MAP by comparing these 
restricted models to the full model (MLM Model 1). MLM Model 2 was constrained so that changes 
in weight and MAP could only be correlated with weight and MAP change in the same and adjacent 
periods, but weight and MAP at baseline were allowed to correlate with changes in weight and MAP in 
all periods. In MLM Model 3 MAP change was not allowed to correlate with weight change in subse- 
quent periods of pregnancy, but weight change was allowed to correlate with later MAP changes, and in 
MLM Model 4 weight change was not allowed to correlate with MAP change in subsequent periods of 
pregnancy, but MAP change was allowed to correlate with later weight changes. These models (MLM 
Models 3 and 4) were used to test whether changes in weight precede changes in MAP or vice versa. 
In MLM Model 5, MAP change up to 18 weeks was not allowed to correlate with weight at baseline or 
weight changes in any period, to test whether baseline weight and changes in weight were only asso- 
ciated with increases in MAP after 18 weeks, but not associated with MAP change prior to 18 weeks, 
when it is decreasing and less likely to be influenced by HDP. Table III shows the deviances of each of 
the models tested. The fits of the constrained models were all much poorer than the full model (MLM 
Model 1) so this was selected as the model that best represented the relationships in the data and all of 
the alternative hypotheses were rejected in favour of this unrestricted model. 

We also fitted an SEM to the data as in Figure 1 with each time point representing a week of preg- 
nancy up to 44 weeks. The earliest week of gestation for which we had data was week 2, producing 43 
time points. We included four splines of gestational age for MAP with the same knot points as in the 
MLM and adjusted for the same set of confounders. All intercept and slope latent variable variances 
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Table III. Model fit comparisons for full and restricted multilevel models adjusted for confounders 
(N = 9429). 


Model 


Deviance 


Difference in df 
(compared with model 1) y^p-value 


MLM Model 1: Full model 


984211.981 




MLM Model 2: Changes in MAP or 
weight are only associated in adjacent 
periods of gestation 


984277.574 


8 < 0.001 


MLM Model 3: Weight change precedes 
MAP change 


984241.507 


3 < 0.001 


MLM Model 4: MAP change precedes 
weight change 


984235.532 


5 < 0.001 


MLM Model 5: Weight at baseline and 
changes in weight are only associated 
with MAP changes after 1 8 weeks 


984228.987 


4 0.002 



Note: 

MLM Model 1: Full model with no restrictions. 



MLM Model 2: Individual-level random effects for periods of MAP and weight change are correlated only with MAP and weight 
changes in concurrent and immediately subsequent periods, but not with later periods (although MAP and weight at 8 weeks (intercept) 
may correlate with MAP and weight changes in any time period). 

MLM Model 3: Individual-level random effects for periods of MAP change are not correlated with weight change in later periods of 
gestation; weight change can correlate with later MAP changes. 

MLM Model 4: Individual-level random effects for periods of weight change are not correlated with MAP change in later periods of 
gestation; MAP change can correlate with later weight changes. 

MLM Model 5: Individual-level random effects for weight at 8 weeks and periods of weight change are not correlated with MAP 
change between 8 and 18 weeks (but may correlate with later MAP changes). 



and covariances were freely estimated, but the residual variances for weight in each time period were 
constrained to be equal, the residual variances for MAP in each time period were constrained to be equal 
and a common residual covariance between weight and MAP at each time point was also estimated, to 
make equivalence with MLM Model 1. Table IV shows the covariances and correlations between the 
intercept and slope latent variables for each of the splines for weight and MAP estimated from the SEM 
in Mplus. The covariances and correlations are very similar to the full MLM Model (MLM Model 1). 
We also compared the estimates of mean weight and MAP at 8 weeks' gestation (baseline) and the mean 
weight and MAP slopes in each period of gestation between multilevel and structural equation models 
that were not adjusted for confounders and these were equivalent in the MLM and SEM. 

Table V shows the associations of weight at 8 weeks and GWG in each period of pregnancy with 
MAP changes in the same and subsequent periods of gestation, with and without adjustment for weight 
and MAP at baseline and earlier changes in these variables. These regression coefficients were derived 
using the variances and covariances of random effects from MLM Model 1, as set out in Equation (7). 
The three different methods for producing the standard errors described in Section 3.4 are compared 
in Table V. We did not use the delta method (Method 2) for regression models that adjusted for more 
than three variables because this method involved specifying an explicit formula for each of the regres- 
sion coefficients and these quickly became extremely complex when adjusting for more variables, and 
impractical to implement. Confidence intervals using Method 1, which assumed the regression coeffi- 
cients had been obtained from a sample of observations of the individual-level random effects, were 
much narrower than for Methods 2 and 3, which treated the regression coefficients as nonlinear trans- 
formations of the estimates of variances and covariances of the random effects, and these two methods 
produced very similar confidence intervals. 

Table V shows that weight at 8 weeks' gestation was positively associated with MAP at 8 weeks and 
with MAP change between 8 and 18 weeks, but negatively associated with MAP change between 18 
and 29 weeks. GWG between 8 and 18 weeks was not strongly associated with MAP in the same period 
but was positively associated with MAP change between 18 and 29 weeks and negatively associated 
with MAP change between 29 and 36 weeks, and these associations remained similar when adjusting 
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for weight and MAP at baseline. GWG between 1 8 and 29 weeks was positively associated with MAP 
changes in the same period but not in subsequent periods and GWG from 29 weeks onwards was posi- 
tively associated with MAP changes between 29 and 36 weeks and from 36 weeks onwards, in models 
with and without adjustment for weight and MAP at baseline and earlier changes in these variables. 

We also considered associations of MAP at 8 weeks and MAP changes with subsequent changes in 
weight, to examine whether there was evidence that MAP increases may influence greater weight gain 
for example through increased oedema (Table SII). In models that were unadjusted for baseline weight 
and earlier GWG and MAP changes, MAP at 8 weeks was negatively associated with GWG in all three 
periods of gestation and MAP changes between 8 and 18 weeks, and between 18 and 29 weeks were 
positively associated with GWG in all subsequent periods. However, after adjustment for weight at 8 
weeks, there was little evidence that MAP at 8 weeks was associated with changes in weight across 
pregnancy, except for a weak association with GWG between 18 and 29 weeks. MAP change between 
8 and 18 weeks was not associated with GWG between 18 and 29 weeks but was weakly positively 
associated with GWG from 29 weeks onwards in models adjusted for weight at 8 weeks and GWG and 
MAP changes in earlier periods of gestation. MAP change between 18 and 29 weeks was not associated 
with subsequent GWG in fully adjusted models. 



6. Discussion 

We have shown that linear splines for two or more change processes can be used in a multivariate 
response MLM to investigate associations between rates of change in adjacent periods and to determine 
whether change in one time period in one variable is associated with change in another variable in the 
next time period, with time periods defined by the splines. It is also possible to adjust for the value of 
one or both of the variables at baseline or changes in the variables in earlier time periods. 

We compared several methods of deriving standard errors for the regression coefficients for the asso- 
ciations between changes in different time periods and found that standard errors using Method 1, 
which acted as though we had a sample of individual-level random effects rather than estimates of their 
variances and covariances, produced confidence intervals that were too narrow and overestimated the 
evidence against the null hypothesis. The delta method is an appropriate way of estimating the variance 
of linear and nonlinear transformations of parameter estimates [27], and incorporated the variance of 
the estimates of the variances and covariances of random effects in the standard errors of the regression 
coefficients derived from these. Using percentiles of the distribution of the regression coefficients to 
combine estimates over 10,000 simulated datasets produced comparable confidence intervals for these 
regression coefficients to those produced using the delta method. However, it should be noted that these 
methods (Methods 2 and 3) both assume that the distributions of the variances and covariances of the 
random effects have a multivariate normal distribution, and so should only be used for large samples 
where this approximation is justified. 

We presented regression coefficients for the associations of pregnancy weight and MAP at baseline 
and rates of change in these variables with subsequent rates of change in the other variable, both with 
and without adjustment for weight and MAP at baseline and changes in earlier periods of pregnancy for 
comparison. This was to assess whether associations between pregnancy weight and MAP change were 
confounded by earlier changes or baseline values. In general, careful consideration is needed to identify 
which variables to adjust for to tease out evidence for causal relationships and this decision should be 
motivated by a particular hypothesis to be tested or by prior knowledge of the relationships between the 
variables studied. 

Previous studies have used covariation in the random effects of two variables in MLMs, either mea- 
sured at a single point in time [35] or longitudinally by examining the covariance between random 
intercept and random slope terms in linear growth models [7, 8] to investigate relationships between two 
variables or rates of change in these variables measured in the same individuals. However, the advantage 
of using linear splines is that the timings of associations can be identified. For this example analysis 
we were able to assess whether changes in weight or MAP preceded changes in the other variable, and 
whether this occurred in all periods of pregnancy or only in early or late pregnancy. This is not possible 
with a bivariate linear growth model. The linear splines also have the advantage of representing nonlin- 
ear patterns of change in an interpretable way, so that the correlations between periods of linear change 
are more meaningful than correlations between random effects for squared or cubic terms in polynomial 
models would be. We could also extend the model to include categorical outcomes, provided that the 
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assumption of normally distributed random effects was not violated by including these variables. It is 
also possible to include different sets of splines for each of the variables if patterns of change in the 
variables are different. 

Cross-lagged SEMs have been used effectively with panel data to demonstrate possible causal relation- 
ships between two or more variables measured at several points in time [36,37] and Granger investigated 
causal relationships between two time series by assessing whether earlier values of each of the vari- 
ables predicted later observations of the other variable [38]. The synchronous common factor model 
may be used when a latent variable is thought to influence the two measured processes simultane- 
ously [39]. These models have the flexibility to include prior information by adding constraints and 
removing unlikely paths from the models. However, the cross-lagged panel design relies on observations 
being made at the same point in time for each individual and for each measure, and these models are 
restricted to researcher-determined timing of measurements, rather than time points that may be more 
biologically meaningful. 

Linear MLMs can be equivalently estimated using the SEM framework with both balanced and unbal- 
anced data structures [31-33] but many nonlinear and generalised linear MLMs cannot be equivalently 
estimated as SEMs [7, 33]. The MLM approach is a more efficient method of fitting linear models with 
unbalanced data because within the MLM framework it is not necessary for data to be measured for all 
individuals at all time points, provided that missingness is at random, whereas in SEMs missing data 
must be dealt with in some way [31] such as by full information maximum likelihood estimation. Also, 
the data management and code become complicated in an SEM when individuals have different numbers 
of measurements and these occur at different times for each individual, and the models are computation- 
ally intensive if there are many possible time points. Furthermore, it is possible to fit a wider range 
of hierarchical models within the MLM framework [31]. However, a particular advantage to the SEM 
framework is its high flexibility to specify complex covariance structures, which may not fit into stan- 
dard MLMs, and the ability to include latent outcome or predictor variables, which are measured by a set 
of indicator variables [32]. It is also possible to investigate mediation by estimating direct and indirect 
effects within the SEM framework and obtain tests of overall model fit, which are not available in MLM 
software [33]. Curran suggested that there was no benefit to using an SEM instead of a standard MLM, 
but that SEMs could be used to incorporate additional model features that could not be included in an 
MLM [32]. The multivariate spline MLMs described here could be extended within the SEM framework 
to produce growth mixture models from the multivariate splines, where patterns of change in two or 
more variables are used as indicators of latent class membership to extract subgroups of the population 
with similar change trajectories. 

The MLMs and maximum likelihood procedure used in the SEM models to deal with missing data 
both assume that the data are missing at random, meaning that the missingness in the data can be 
explained by observed variables [40]. If there is reason to suspect that the missing data mechanism 
is nonignorable, there are methods available to test the sensitivity of inferences to data being not missing 
at random for multivariate longitudinal data [41]. 

We have demonstrated that these models can be fitted in both MLM and SEM frameworks. The SEM 
took longer to converge than the MLM, because of the low coverage of the data. All of the women had 
far fewer weight and MAP observations than the number of time periods used and the SEM treated each 
time period in which the individual did not have a measurement as a missing value that needed to be 
estimated using maximum likelihood estimation. However, despite the difficulties with estimation in the 
SEM, the variance-covariance matrix of the latent intercept and slope factors produced by this model 
was very similar to that of the individual-level random effects produced by the MLM. 

Convergence can also be an issue in MLMs, especially if there are many splines and therefore many 
variances and covariances of individual level random effects to be estimated. If the model does not 
converge or calculates negative variances for some of the random effects, it is not possible to derive 
regression coefficients from these. The choice of the number of splines to use for each variable should 
depend on the pattern of change in the variable over time, and should be sufficient to achieve an accept- 
able fit to the observations. It is possible to include a greater number of splines if there is a larger sample 
size and a high number of measurements per individual, and convergence may be easier if the vari- 
ables are measured at the same time points for each individual and fewer covariates are included in the 
model. Difficulties in achieving convergence could be overcome by adding one random effect at a time 
to the model and reestimating using initial values from the previous model, using fewer splines, or by 
reducing the accuracy of the timings of the observations, for example rounding gestational age to whole 
weeks rather than the nearest day, because this reduces the variability in measurement times between 
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individuals. The complexity of the model fitted in our example appeared to be close to the limit of what 
is possible with this dataset, as in some of our constrained models convergence became difficult. 

In summary, we found that GWG up to 18 weeks' gestation was associated with a greater increase in 
MAP in the next period of pregnancy (18 — 29 weeks), but associated with a smaller rise in MAP later 
on (29 — 36 weeks). GWG between 18 and 29 weeks was positively associated with concurrent increases 
in MAP but not with MAP in later periods, and GWG from 29 weeks onwards was also associated with 
concurrent MAP. There was some evidence that MAP change was also related to later GWG, with a 
weak positive association between MAP change up to 18 weeks and GWG from 29 weeks onwards. 
This method could be applied to many other areas where it is of interest whether changes in one variable 
precede changes in another, and in what time periods these relationships exist and are strongest. For 
example it could be used to investigate whether changes in weight are related to changes in respiratory 
function or bone mass or whether changes in socioeconomic position over time precede changes in a 
wide range of health and social outcomes. 
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